home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
SGI Hot Mix 17
/
Hot Mix 17.iso
/
HM17_SGI
/
research
/
lib
/
linfit.pro
< prev
next >
Wrap
Text File
|
1997-07-08
|
6KB
|
150 lines
;$Id: linfit.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
;
; Copyright (c) 1994-1997, Research Systems, Inc. All rights reserved.
; Unauthorized reproduction prohibited.
;+
; NAME:
; LINFIT
;
; PURPOSE:
; This function fits the paired data {X(i), Y(i)} to the linear model,
; y = A + Bx, by minimizing the chi-square error statistic. The result
; is a two-element vector containing the model parameters,[A,B].
;
; CATEGORY:
; Statistics.
;
; CALLING SEQUENCE:
; Result = LINFIT(X, Y)
;
; INPUTS:
; X: An n-element vector of type integer, float or double.
;
; Y: An n-element vector of type integer, float or double.
;
; KEYWORD PARAMETERS:
; CHISQ: Use this keyword to specify a named variable which returns the
; chi-square error statistic as the sum of squared errors between
; Y(i) and A + BX(i). If individual standard deviations are
; supplied, then the chi-square error statistic is computed as
; the sum of squared errors divided by the standard deviations.
;
; DOUBLE: If set to a non-zero value, computations are done in double
; precision arithmetic.
;
; PROB: Use this keyword to specify a named variable which returns the
; probability that the computed fit would have a value of CHISQR
; or greater. If PROB is greater than 0.1, the model parameters
; are "believable". If PROB is less than 0.1, the accuracy of the
; model parameters is questionable.
;
; SDEV: An n-element vector of type integer, float or double that
; specifies the individual standard deviations for {X(i), Y(i)}.
;
; SIGMA: Use this keyword to specify a named variable which returns a
; two-element vector of probable uncertainties for the model par-
; ameters, [SIG_A,SIG_B].
;
; EXAMPLE:
; Define two n-element vectors of paired data.
; x = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89, -0.12, 1.41, $
; 2.95, 2.18, 3.72, 5.26]
; y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98, -2.87, -1.66, $
; -0.78, -2.61, 0.31, 1.74]
; Define a vector of standard deviations with a constant value of 0.85
; sdev = replicate(0.85, n_elements(x))
; Compute the model parameters, A and B.
; result = linfit(x, y, chisq = chisq, prob = prob, sdev = sdev)
; The result should be the two-element vector:
; [-3.44596, 0.867329]
; The keyword parameters should be returned as:
; chisq = 11.4998, prob = 0.319925
;
; REFERENCE:
; Numerical Recipes, The Art of Scientific Computing (Second Edition)
; Cambridge University Press
; ISBN 0-521-43108-5
;
; MODIFICATION HISTORY:
; Written by: GGS, RSI, September 1994
; LINFIT is based on the routines: fit.c, gammq.c, gser.c,
; and gcf.c described in section 15.2 of Numerical Recipes,
; The Art of Scientific Computing (Second Edition), and is
; used by permission.
; Modified: SVP, RSI, June 1996
; Changed SIG_AB to SIGMA to be consistant with the other
; fitting functions. Changed CHISQR to CHISQ in the docs
; for the same reason. Note that the chisqr and the SIG_AB
; keywords are left for backwards compatibility.
; Modified: GGS, RSI, October 1996
; Modified keyword checking and use of double precision.
; Added DOUBLE keyword.
;-
FUNCTION LinFit, x, y, chisqr = chisqr, Double = Double, prob = prob, $
sdev = sdev, sig_ab = sig_ab, sigma = sigma
ON_ERROR, 2
TypeX = SIZE(X)
TypeY = SIZE(Y)
nX = TypeX[TypeX[0]+2]
nY = TypeY[TypeY[0]+2]
if nX ne nY then $
MESSAGE, "X and Y must be vectors of equal length."
;If the DOUBLE keyword is not set then the internal precision and
;result are identical to the type of input.
if N_ELEMENTS(Double) eq 0 then $
Double = (TypeX[TypeX[0]+1] eq 5 or TypeY[TypeY[0]+1] eq 5)
nsdev = n_elements(sdev)
if nsdev eq nX then begin ;Standard deviations are supplied.
wt = 1.0 / sdev^2
ss = TOTAL(wt, Double = Double)
sx = TOTAL(wt * x, Double = Double)
sy = TOTAL(wt * y, Double = Double)
t = (x - sx/ss) / sdev
st2 = TOTAL(t^2, Double = Double)
b = TOTAL(t * y / sdev, Double = Double)
endif else if nsdev eq 0 then begin
ss = nX + 0.0
sx = TOTAL(x, Double = Double)
sy = TOTAL(y, Double = Double)
t = x - sx/ss
st2 = TOTAL(t^2, Double = Double)
b = TOTAL(t * y, Double = Double)
endif else $
MESSAGE, "sdev and x must be vectors of equal length."
if Double eq 0 then begin
st2 = FLOAT(st2) & b = FLOAT(b)
endif
b = b / st2
a = (sy - sx * b) / ss
sdeva = SQRT((1.0 + sx * sx / (ss * st2)) / ss)
sdevb = SQRT(1.0 / st2)
if nsdev ne 0 then begin
chisqr = TOTAL( ((y - a - b * x) / sdev)^2, Double = Double )
if Double eq 0 then chisqr = FLOAT(chisqr)
prob = 1 - IGAMMA(0.5*(nX-2), 0.5*chisqr)
endif else begin
chisqr = TOTAL( (y - a - b * x)^2, Double = Double )
if Double eq 0 then chisqr = FLOAT(chisqr)
prob = chisqr * 0 + 1 ;Make prob same type as chisqr.
sdevdat = SQRT(chisqr / (nX-2))
sdeva = sdeva * sdevdat
sdevb = sdevb * sdevdat
endelse
sig_ab = [sdeva, sdevb]
sigma = sig_ab
if Double eq 0 then RETURN, FLOAT([a, b]) else RETURN, [a, b]
END